# AMAR ET AL. - COUNTERING MISINFORMATION EARLY (2025)
## REPLICATION FILE: 04_compliance.R
### This script models compliance predictors and interactions with treatment assignment.
# ----
# Compliance predictors ----
# set relevant predictors
compliance_predictors <- list(
  list(iv = "gender", label = "Gender", type = "Demographic"),
  list(iv = "class", label = "Grade", type = "Demographic"),
  list(iv = "age", label = "Age", type = "Demographic"),
  list(iv = "religion_hindu_num", label = "Religion - Hindu", type = "Demographic"),
  list(iv = "language_hindu_num", label = "Language - Hindi", type = "Demographic"),
  list(iv = "asset_index", label = "Asset Index", type = "Demographic"),
  list(iv = "fathers_education", label = "Father's Education", type = "Education"),
  list(iv = "mothers_education", label = "Mother's Education", type = "Education"),
  list(iv = "school_gov_num", label = "Government School", type = "Education"),
  list(iv = "science", label = "Science Knowledge", type = "Education"),
  list(iv = "mobile_internet_num", label = "Mobile Internet", type = "Demographic"),
  list(iv = "trust_newspapers_num", label = "Newspapers", type = "Trust"),
  list(iv = "trust_social_media_num", label = "Social Media", type = "Trust"),
  list(iv = "trust_tv_num", label = "TV", type = "Trust"),
  list(iv = "trust_friends_family_num", label = "Friends and Family", type = "Trust"),
  list(iv = "vaccinated_num", label = "Vaccinated", type = "Trust"),
  list(iv = "ayurveda_effective_num", label = "Ayurveda", type = "Trust"))

compliance_predictors <- rbindlist(compliance_predictors, fill = TRUE, use.names = TRUE)

compliance_predictors <- data.frame(compliance_predictors)

# function to tabulate regressions for each dv, with baseline controls
tab_regression_compliance <- function(dv, iv, design) {
  fml <- as.formula(paste(dv, "~", iv, "+ library_spillover_pre"))  
  mod <- svyglm(fml, design = design)
  out <- tidy(mod) %>%
    mutate(
      dv = dv, 
      iv = iv,
      n = nobs(mod),
      .before = term
    )
  out
}


regressions_compliance <- lapply(compliance_predictors$iv, function(iv) {
  tab_regression_compliance("compliance_minimal", iv, bimli_svy_dkr.rm)
}) %>%
  bind_rows() %>%
  mutate(
    margin = qnorm(0.975) * std.error,
    lower = estimate - margin,
    upper = estimate + margin
  ) %>%
  filter(term != "(Intercept)" & !grepl("library_spillover_pre", term)) # remove intercepts and district FEs

regressions_compliance <- left_join(compliance_predictors, 
                                    regressions_compliance, by = "iv") %>%
  rowwise() %>%
  mutate(
    label_rec = case_when(grepl(iv, term) & 
                            nchar(term) > nchar(iv) ~ 
                            paste(label, "-", sub(iv, "", term)),
                          TRUE ~ label)) %>%
  ungroup()

## Plot ----
regressions_compliance %>%
  filter(iv != "treatment") %>%
  mutate(significant = if_else(p.value <= 0.05, "Yes", "No")) %>%
  ggplot(aes(x = estimate, y = fct_rev(label_rec), color = significant)) +
  geom_vline(xintercept = 0, lty = 2) +
  geom_point() +
  geom_errorbar(aes(xmin = lower, xmax = upper), width = 0) +
  facet_grid(rows = "type", scales = "free_y", space = "free_y", switch = "y") +
  xlab("Estimated effect on compliance (# of sessions)") +
  scale_y_discrete(position = "right") +
  scale_color_manual(values = c("black", "red")) +
  theme_bw() %+replace%
  theme(axis.title.y = element_blank(),
        strip.text.y.left = element_text(angle = 0),
        legend.position = "none")
ggsave("output/figures/compliance_predictors.pdf", height = 5, width = 8)

## Table ----
regressions_compliance %>%
  #filter(idx == type & secondary_outcome == FALSE) %>% # filter to variables that are included in the index
  select(label_rec, n, estimate, std.error, p.value) %>%
  mutate(n = format(n, big.mark = ","),
         sig.stars = case_when(p.value < 0.001 ~ "***",
                               p.value < 0.01 ~ "**",
                               p.value < 0.05 ~ "*",
                               TRUE ~ ""),
         estimate = str_c(format(round(estimate, 2), nsmall = 2), sig.stars), # add significance stars to estimate
         std.error = round(std.error, 3),
         p.value = case_when(p.value > 0.9 ~ "$>$0.9",
                             p.value < 0.001 ~ "$<$0.001",
                             p.value < 0.01 ~ as.character(round(p.value, 3)),
                             TRUE ~ as.character(round(p.value, 2)))) %>%
  select(-sig.stars) %>%
  rename(Predictor = label_rec,
         N = n,
         Estimate = estimate,
         SE = std.error,
         "p-value" = p.value) %>%
  xtable::xtable(caption = "Compliance predictors", align = "llllll", digits = 3,
                 label = "tab:compliance") %>%
  print(include.rownames = FALSE, 
        caption.placement = "top",
        sanitize.text.function = identity, # add backslashes to special characters
        sanitize.colnames.function = \(x) str_c("\\textbf{", x, "}"),
        add.to.row = list(pos = list(nrow(.)), 
                          command = paste(
                            "\\hline\n\\multicolumn{5}{l}{\\footnotesize *p$<$0.05; **p$<$0.01; ***p$<$0.001. Models include library-spillover strata FEs.} \\\\\n")),
        hline.after = c(-1,0), # add horizontal lane after second to last column
        comment = FALSE,
        table.placement = getOption("xtable.table.placement", "h!"),
        file = paste0("output/tables/tab_compliance", ".tex"))

# Compliance predictors by condition ----
# function to tabulate regressions for each dv, with baseline controls
tab_regression_compliance_int <- function(dv, iv, design) {
  fml <- as.formula(paste(dv, "~", iv, "+", "treatment+", iv, "* treatment", "+ library_spillover_pre"))  
  mod <- svyglm(fml, design = design)
  out <- tidy(mod) %>%
    mutate(
      dv = dv, 
      iv = iv,
      n = nobs(mod),
      .before = term
    )
  out
}

compliance_predictors_int <- compliance_predictors %>%
  filter(iv != "treatment") 

regressions_compliance_int <- lapply(compliance_predictors_int$iv, function(iv) {
  tab_regression_compliance_int("compliance_minimal", iv, bimli_svy_dkr.rm)
}) %>%
  bind_rows() %>%
  mutate(
    margin = qnorm(0.975) * std.error,
    lower = estimate - margin,
    upper = estimate + margin
  ) %>%
  filter(term != "(Intercept)" & !grepl("library_spillover_pre", term) &
           grepl(":treatment", term))


regressions_compliance_int <- left_join(compliance_predictors_int, regressions_compliance_int, by = "iv") %>%
  rowwise() %>%
  mutate(
    label_rec = case_when(grepl(iv, term) & nchar(term) > nchar(iv) + nchar(":treatmentMedia Literacy") ~ paste(label, "-", sub(iv, "", term)),
                          TRUE ~ label),
    label_rec = sub(":treatmentMedia Literacy", "",label_rec),
    label_rec = paste0(label_rec, " * T")) %>%
  ungroup()

## Plot ----
regressions_compliance_int %>%
  mutate(significant = if_else(p.value <= 0.05, "Yes", "No")) %>%
  ggplot(aes(x = estimate, y = fct_rev(label_rec), color = significant)) +
  geom_vline(xintercept = 0, lty = 2) +
  geom_point() +
  geom_errorbar(aes(xmin = lower, xmax = upper), width = 0) +
  facet_grid(rows = "type", scales = "free_y", space = "free_y", switch = "y") +
  xlab("Interaction effect of covariates and treatment on compliance") +
  scale_y_discrete(position = "right") +
  scale_color_manual(values = c("black", "red")) +
  theme_bw() %+replace%
  theme(axis.title.y = element_blank(),
        strip.text.y.left = element_text(angle = 0),
        legend.position = "none")
ggsave("output/figures/compliance_predictors_interaction.pdf", height = 5, width = 8)

## Table ----
regressions_compliance_int %>%
  select(label_rec, n, estimate, std.error, p.value) %>%
  mutate(n = format(n, big.mark = ","),
         sig.stars = case_when(p.value < 0.001 ~ "***",
                               p.value < 0.01 ~ "**",
                               p.value < 0.05 ~ "*",
                               TRUE ~ ""),
         estimate = str_c(format(round(estimate, 2), nsmall = 2), sig.stars), # add significance stars to estimate
         std.error = round(std.error, 3),
         p.value = case_when(p.value > 0.9 ~ "$>$0.9",
                             p.value < 0.001 ~ "$<$0.001",
                             p.value < 0.01 ~ as.character(round(p.value, 3)),
                             TRUE ~ as.character(round(p.value, 2)))) %>%
  select(-sig.stars) %>%
  rename(Predictor = label_rec,
         N = n,
         Estimate = estimate,
         SE = std.error,
         "p-value" = p.value) %>%
  xtable::xtable(caption = "Compliance predictors (*treatment)", align = "llllll", digits = 3,
                 label = "tab:compliance_int") %>%
  print(include.rownames = FALSE, 
        caption.placement = "top",
        sanitize.text.function = identity, # add backslashes to special characters
        sanitize.colnames.function = \(x) str_c("\\textbf{", x, "}"),
        add.to.row = list(pos = list(nrow(.)), 
                          command = paste(
                            "\\hline\n\\multicolumn{5}{l}{\\footnotesize *p$<$0.05; **p$<$0.01; ***p$<$0.001. Models include library-spillover strata FEs.} \\\\\n")),
        hline.after = c(-1,0), # add horizontal lane after second to last column
        comment = FALSE,
        table.placement = getOption("xtable.table.placement", "h!"),
        file = paste0("output/tables/tab_compliance_int", ".tex"))

# END of 04_compliance.R ----
